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Summary 

This paper is devoted to show how to use automatic 
differentiation in reverse mode as a powerful tool in 
optimization procedures. It is also shown that for 
aerodynamic applications the gradients have to be 
as accurate as possible. In particular, the effect of 
having the exact gradient of the first or second or- 
der spatial discretization schemes is presented. We 
show that the loss of precision in the gradient affects 
not only the convergence, but also the final shape. 
Both two and three dimensional configurations of 
transonic and supersonic flows have been investi- 
gated. These cases involve up to several thousand 
of control parameters. 

1 Introduction 

The problems of interest here belong to optimal 
shape design in aeronautics. Some implicit cost 
functional has to be minimized over a set of possible 
states, under the constraint that the state equations 
(steady Euler) are satisfied. 

When using gradient based methods for optimiza- 
tion, we need the gradient of the cost function with 
respect to control parameter variations. This is a 
severe limitation if a direct method based on suc- 
cessive cost function evaluation is used. Indeed, the 
cost of one evaluation of the gradient will be pro- 
portional to the number of control parameters. 

We use the Odyssee system [1,2] in reverse mode 
to produce the Jacobian of the cost function with 
respect to the control variables. Odyssee takes as 
input a Fortran 77 program and a set of variables 
and returns a new Fortran 77 program computing 
the derivatives of the original function with respect 
to the given variables. This gradient has been used 
in a projected conjugate gradient method for mini- 
mization. 

As mentioned before, in optimization procedures 
we have to compute the effects of control variables 


variations on the cost J (i.e. dJ/dx c ). Therefore, 
if a control point is moved, this variation has to 
be propagated to all the mesh nodes. Because the 
target geometries are usually described by an un- 
structured mesh, we have developped a correspond- 
ing framework for shape deformations to take into 
account these deformations over the meshes. The 
meshes used here contain only triangles in 2D and 
tetraedra in 3D. 

2 Control Problem Formula- 
tion 

In this paper the aim is to minimize a functional 
J(x 1 U(x)) under geometric and aerodynamic con- 
straints. Here, x indicates the geometrical de- 
scription of a configuration and the flow pattern 
around this shape is the solution of the steady Eu- 
ler system of fluid dynamics in conservation form: 
V.(F(U)) = 0, where U is the vector of conserva- 
tive variables (i.e. U = (p, pu, p(C v T + ^M 2 ))*), F 
represents the advective operator. This system has 
4 equations in 2D (5 equations in 3D) for 5 variables 
(6 variables in 3D) and the system is closed using 
the equation of state p = p(p, T). 

2.1 The Reverse Mode 

To calculate the gradient of J under aerodynamic 
constraints, automatic differentiation in reverse 
mode has been used. In this approach, the lines 
of the programs describing the relations between 
the variation of the design variables and the cost 
function including the grid and the ’steady’ flow 
equations are multiplied by parameters (p) and an 
augmented Lagrangian (L) is constructed. The val- 
ues of the parameters are obtained from the opti- 
mality conditions (i.e. that the first variations of L 
with respect to intermediate variables vanish). The 
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solution can be always obtained simply by back sub- 
stitution (hence the notion of reverse mode) . Once 
the parameters are evaluated, the gradient of L can 
be easily calculated (a simple example is given in 
the appendix). 

In this approach no adjoint system is solved and 
the work in the back substitution step is equivalent 
to one explicit iteration of the governing equations 
[3], 

The advantages over the other methods are clear. 
In particular, this leads to the exact gradient of the 
discrete cost function with respect to the control 
variables and the computational time for the gradi- 
ent is independant of the number of controls. 

This gradient is then used in a conjugate gradi- 
ent method to solve the optimization problem. In 
[4,5,6], the gradients obtained by Odyssee have been 
compared with those obtained using finite differ- 
ences for similar problems. 

2*2 The Conjugate Gradient Method 

Our minimization tool is quite simple. It is based on 
a conjugate gradient method with optimal descent 
step. We use projection to take into account the 
local geometrical constraints and global constraints 
are taken into account in the cost function using 
Lagrange multipliers. 

More precisely, the algorithm is as follows (we de- 
note J(z,U(x)) by J{z): 

zogiven , 

do 1 iteration of the steepest descent algorithm, 
for n = 2, 3, ... do 

. I 

7 |v r j(*"- 2 )P’ 

h n = -V r J(x n - 1 ) + 7 /» n " 1 , 


2.3 Flow Solver 

The NSC2KE solver uses a finite element/finite vol- 
ume formulation on unstructured meshes involving 
triangles in 2D and tetraedra in 3D. Second or- 
der accuracy in space has been achieved using a 
MUSCL type reconstruction and limiters have been 
used to prevent oscillations. The time dependant 
equation (dU/dt + V.F(U) = 0) is marched in time 
to a steady state. The time discretization is based 
on a 4-stage Runge-Kutta scheme. We will show 
that for the optimization procedure it is important 
to have a gradient including all of these ingredients 
(especially, the second order MUSCL reconstruction 
step). 

2*4 Geometry Modifications and 
Control Points Definition 

In 2D, the control points are fitted by a cubic 
spline. The splines have two features. They have 
a smoothing effect on the variations of the control 
points (< 5z c ) and they propagate these variations to 
the other body points which are not control points 

For 3D applications, the use of generalized surface 
splines is quite complicated and involves CAD con- 
cepts and deriving these objects is more difficult 
than the fluid solver. The present unstructured 
framework [3] for geometry modifications is based 
on the following: 

1. All the nodes on the shape are control points. 

2. To avoid oscillations, a smoothing operator is 
defined over the shape. This can be, for instance, a 
few Jacobi iterations to solve ((I — eA)Sx w = Sx w ), 
where Sz w is the smoothed shape variation for the 
shape nodes and Sz w is the variation given by the 
optimization tool. Once (x w ) known, we have to 
expand these variations overall the mesh. This is 
done by solving an elliptic system. These tools have 
also been derived by the automatic differentiation 
procedure. 


x n = x"- 1 + A"fc n , 

where, 

A n = mm* J(x, t/(x)). 
with x = x n_1 - AV r J(x n_1 ) and 


„ j oj t dJ. T( du. 

V '' = s + <»Z7> <&r>- 


This algorithm converges to a local minimum of J . 


2*5 Geometrical Constraints 

The present geometrical constraints are of two 
types. The first one is imposed by defining two lim- 
iting surfaces (curves in 2D) between which shape 
variations are allowed. As all shapes in this paper 
are of wing type, the second constraint is that the 
original planform should remain unchanged. This 
means for instance in 2D that the leading and trail- 
ing edge are frozen. 
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3 Results 

Two and three dimensional results of inverse prob- 
lems and drag reduction for airfoils and wings are 
presented for transonic and supersonic flows. 


3.1 An Inverse Problem in 2D 

The first case consists of an inverse problem based 
on a given pressure distribution. The cost function 
is given by J(x) = \ f g \p x - p ta rgtt\ 2 dx, where, 
Ptarget is a given target pressure and p T the actual 
flow pressure. The design takes place at Mach num- 
ber 0.85 and zero angle of incidence. The initial 
shape is the RAE2822. The target pressure corre- 
sponds to the same airfoil deformed by about 30 
percent on the upper surface and 20 percent on the 
lower surface. This leads to a shift of the shocks to 
the right by about 20 percent of the chord. There 
are 20 control points on the airfoil and 60 total 
nodes. This is a quite coarse mesh and enables us 
to compare the Odyssee gradients with finite differ- 
ences (see Fig. 1). 

The cost function has been computed using a second 
order scheme but the gradient has been computed 
either by a first or a second order scheme. This is 
to show that this loss of precision for the gradient 
impacts both the convergence and the final shape. 


3.2 A drag reduction problem in 2D 

The aim here is to reduce the shock-induced drag for 
a RAE2822 profile at Mach number 0.8 and zero an- 
gle of incidence. The shock on the upper surface is 
quite close to the trailing edge and is more difficult 
to remove as in this region geometrical constraints 
are more important. On the other hand, the shock 
on the lower surface is easier to remove. 

As in [7,8], consider as cost function J(x) = 
\ Jr |p — p, | 2 dx-h lOCrf, where pi is the original pres- 
sure distribution and Cd the drag coefficient. The 
first term forces the profile to remain as much as 
possible close to the original shape. Moreover, the 
shape variation is limited to 5 percent. The drag 
(resp. lift) has been reduced from 2.110“ 2 (resp. 
0.292) to 1.2510"* 2 (resp. 0.291) with the second 
order gradient and to 1.5310” 2 (resp. 0.289) using 
the first order one. The airfoil volume’s has almost 
remained unchanged. 


3.3 An Inverse Problem in 3D 

We consider a pressure recovery problem in 3D. The 
original shape here is the ONERA M6 wing. The 
target shape is the same wing deformed on its upper 
surface by 30 percent. Our aim is to recover this 
shape using the corresponding pressure distribution 
starting from the M6 wing. The section definitions 
of the wing are not used. Also, the cross sections 
have been computed by interpolation (see Fig. 6). 
All the wall mesh points are control points. In this 
case, we have about 700 nodes on the wing. It is 
therefore impossible to treat this case without the 
adjoint (inverse) mode. The important remark here 
is that when we use only first order gradient, the 
target shape is not correctly recovered. 

3.4 Wave-Drag Reduction in 3D 

The sum here is to reduce the shock-induced drag 
over an M6 type wing. The original wing is an M6 
wing with the upper surface deformed by 10 per- 
cent to obtain a nonsymmetric wing. The mesh 
has around 10 5 tetraedra and there are about 2000 
control points on the wing (see. Fig. ii). The 
design take place at Mach number 0.84 and in- 
cidence of 3.06. This configuration involves a A- 
shock on the upper surface and our aim is to pro- 
duce a wing as close as possible to the original with 
smoothed shocks. We will use the second order 
gradient in this computation. The cost function is 
J(x) = \ Jr IF ” Po\ 2 dx + 10 Cd + | Ci - Cf|, where 
Ci and C? are the actual and initial lift coefficients. 
The drag has been reduced by about 10 percent. 

3.5 A supersonic case 

This is an inverse problem at Mach 3 over an 
Naca0014 airfoil. The aim is to show that the ex- 
tension of our approach to supersonic flows does not 
introduce any particular difficulty even at bound- 
aries. 


4 Concluding Remarks 

A new approach involving the reverse mode of auto- 
matic differentiation has been presented. A general 
framework for treating geometries with unstruc- 
tured discretization has been introduced for both 
two and three dimensional configurations. Prelim- 
inary examples show the ability of the method to 
treat inverse and control problems. 
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The inverse mode is shown to be a powerful tool 
for providing the exact gradient of the discrete op- 
erators. Thus, 3D cases with several thousands of 
control points are possible to calculate. 

It has also been pointed out that the gradient 
should include all the ingredients of the discrete op- 
erators. The accuracy of the gradient impacts not 
only the convergence but also the final shapes. 

It is important to notice that, except for the gra- 
dient which use a backward time integration pro- 
cedure, all the calculations have been done by ex- 
plicit schemes. A parallel implementation of this 
approach is therefore quite simple. 

Future works will include more extensive validation 
of these techniques in 3D configurations. 
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and the design variables x such that f(u(x)) y where 
ti(x) is the solution of the flow equations. Assume 
the following Fortran 77program: 

ul=x 

u2=x**2+2*ul 

f=ul+u2 

In automatic differentiation in reverse mode, we 
consider the lines of the program as constraints and 
associate to each of them a Lagrange multiplier and 
define an augmented Lagrangian as follows: 

L = til + «2 + Pl(«2 ” z 2 - 2lli) -f p 2 (u I - x). 


We know that at the solution we have: 


- — = 1 — 2pi + p 2 = 0, 
OU 1 


8L 

du 2 


1 +Pi = 0. 


We notice that to find p,-, we have to solve the pre- 
vious set of equations in ’’reverse” order. Once pi 
known, we have: 


dL_ 

dx 


5/ 

dx 


-2 Pix - P2, 


which is the Jacobian of /. 



Figure 1: Inverse problem with inverse mode: com- 
parison of the gradients obtained by Odyssee and 
finite differences. 


6 Appendix 

We give a simple example of the automatic differen- 
tiation in reverse mode. Consider a cost function / 
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Figure 2: Inverse problem with inverse mode: ini- 
tial, target and computed shapes using the first and 
second order Jacobians. 


COSTf WU/TDN MBTCWY 



Figure 3: Inverse problem with inverse mode: con- 
vergence histories for the cost using first and second 
order gradients. 



Figure 4: Drag reduction: convergence histories for 
the cost function using second order operator and 
first or second order gradients . 



Figure 5: Drag reduction: initial and final shapes 
obtained with the first and second order gradients. 



Figure 6: 3D inverse problem: M6 wing, upper sur- 
face discretization. 



Figure 7: 3D inverse problem: Iso-Mach contours 
over the M6 wing } s upper surface (in the range 
[0.3, 1.4] with 40 contours). 
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Figure 8: 3D inverse problem : Iso-Mach contours 
over the target wing ’ s upper surface (in the range 
[0,3, 1.4] 40 contours ). 



<1 » U M U dl I It 


Figure 9: 3D inverse problem: Iso-Mach contours 
over the upper surface of the shape obtained using 
the second order gradient (uniform discretization of 
the interval [0.3, 1.4] with 40 isovalues). 


Figure 11: 3D Drag reduction: nonsymmetrxc wing, 
upper surface discretization , all these points are 
control points. 



Figure 12: 3D drag reduction: Iso-Mach contours 
over the initial wing’s upper surface. 



Figure 10: 3D inverse problem: Iso-Mach contours 
over the upper surface of the shape obtained using 
the first order gradient (in the range [0.3,1. 4] with 
40 contours). 
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Figure 13: 3D drag reduction: Iso-Mach contours 
over the upper surface of the shape obtained using 
the second order gradient. 











Figure 14: 3D drag reduction: shape cross-section 
for the initial and the optimized shape after 5 design 
iterations. 




Figure 17: Supersonic flow at M= 3. Pressure dis- 
tribution: initial , target and computed. 


Figure 15: 3D drag reduction: pressure cross- 

section over the initial and the optimized shape after 
5 design iterations. 




Figure 18: Supersonic flow at M~3. Initial , target 
and final shapes. 


Figure 16: Supersonic flow at M=3. Convergence 
history. 
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